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ABSTRACT 

We investigate the formation and properties of galaxies hosting z ~ 6 quasars, in 
the gigaparsec scale cosmological hydro dynamical simulation: MassiveBlack , which 
includes a self-consistent model for star formation, black hole accretion and associated 
feedback. We show that the MassiveBlack reproduces current estimates of the galaxy 
stellar mass function z = 5,6. We find that quasar hosts in the simulation are compact 
gas rich systems with high star formations rates of SFR ~ 100 — lO 3 M yr _1 consistent 
with observed properties of Sloan quasar hosts in the redshift range 5.5 < z < 6.5. 
We show that the star-forming gas in these galaxies predominantly originates from 
high density cold streams which are able to penetrate the halo and grow the galaxy at 
the center. MassiveBlack predicts a deviation from the local M BH — a and M BH — M* 
relation implying that black holes are relatively more massive for a given stellar host 
at these redshifts. 

Key words: galaxies: active - galaxies: formation - galaxies: evolution - methods: 
numerical quasars: general - black hole physics 



1 INTRODUCTION 

Supermassive black holes (SMBH) are now ubiquitously 
found in the nuclei of local galaxies. Tight correlations have 
been observed between the central black hole and its host 
galaxy (Magorrian et al. 1998; Ferrarese & Merritt 2000; 
Gebhardt et al. 2000; Graham et al. 2001; Tremaine et al. 
2002; Haring & Rix 2004; Marconi & Hunt 2003), imply- 
ing that the growth of the quasar, powered by the SMBH, 
is intimately linked to the formation of its host galaxy. At 
z ~ 6 and above the most direct constraint on the evolution 
of SMBHs comes from observations of luminous quasars in 
the Sloan Digital Sky Survey (Fan et al. 2006; Jiang et al. 
2009). Recently a detection has been confirmed for a quasar 
at z = 7 (Mortlock et al. 2011). These quasars are mainly 
optically selected and represent the bright end of the quasar 
population at this epoch. They are rare with number den- 
sities of n rsj a few Gpc -3 and extremely luminous with 
inferred masses of M BH ~ 1O 9 M , suggesting that they are 
harboured in rare halos of mass Mhaio ~ 1O 13 M at these 
redshifts. Observations in other bands, far infrared (FIR)- 
radio, suggest that the spectral energy distributions (SED) 
of these quasars are similar to those of their low redshift 
counterparts (Wang et al. 2008). This means that 1O 9 M 
black holes are fully developed and already in place even 
when the Universe was relatively young (< 10 9 years at 
z > 6). 

Reprocessed thermal dust continuum emission in the 
FIR provides a clean method for deriving the total star for- 



mation rates in galaxies (Dale & Helou 2002). Bright quasars 
also add to the FIR luminosity, Lfir, and one needs to cor- 
rect for it to estimate the star formation rate (SFR) of the 
host galaxy. The excess Lfir (corrected) for a sample of 
z rsj 6 quasar hosts suggests a massive starburst origin (SFR 
- 10 2 - 7 - 10 3 4 M o yr _1 ) for them (Wang et al. 2010, 2011). 
These authors also detected large reservoirs of molecular 
gas, with mass M mol ~ 1O 1O M , in these galaxies through 
the emission of redshifted carbon monoxide (CO) line, fur- 
ther corroborating the result that the observed large star 
formation activity in these galaxies are sustained through 
an abundant supply of molecular gas, the fuel for star for- 
mation. The cold gas in these galaxies is localised within 
a radius of a few kpc (Walter et al. 2004, 2009; Wang et 
al. 2010). The relation between CO luminosity, £00(1-0)5 
and Lfir for these galaxies are similar to those of typical 
star-forming systems at lower redshift. 

However the relation between the black hole mass, M BH 
and the bulge velocity dispersion of its host, a, is seen to 
be above the local relation(Wang et al. 2010), indicating 
that the M BH — a relation is evolving. This is true even 
when some of the assumptions, like degeneracy of a with 
inclination angle are considered. Even though there is still 
significant debate whether there may be observational biases 
influencing these results (Lauer et al. 2007). 

Recent observations with the refurbished Hubble Space 
Telescope (HST) have yielded a considerably larger sam- 
ple of more typical galaxies at 5 < z < 8. The newly in- 
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stalled Wide Field Camera (WFC3/IR) and the Advanced 
Camera for Surveys (ACS) on the HST have been used to 
detect much fainter galaxies out to z ~ 8. These provide 
strong constraints on the UV luminosity function of galax- 
ies (Bouwens et al. 2011a,b) at these redshifts. In addition 
WFC3/IR+ACS in combination with the Spitzer telescope 
has been used to detect the highest redshift (z ~ 10) galaxy 
till date (Bouwens et al. 2011c; Oesch et al. 2011). The UV 
luminosity function at these redshifts has been used to esti- 
mate the galaxy stellar mass function (GSMF) from z ~ 5—7 
(Gonzalez et al. 2011). 

These observations of high redshift galaxies impose 
strong constraints on theoretical models of galaxy formation. 
Hydrodynamic simulations (with models for star formation) 
have been carried out for studying the global properties of 
high redshift galaxies (Salvaterra, Ferrara, & Dayal 2011; 
Jaacks, Choi, h Nagamine 2011). They are able to reproduce 
the observed UV luminosity of these galaxies with some as- 
sumptions for extinction due to dust. However the shape of 
the GSMF (Jaacks, Choi, & Nagamine 2011) is found to be 
inconsistent with those inferred from observations (Gonzalez 
et al. 2011). These simulations do not model the growth of 
black holes and their feedback on the surrounding environ- 
ment and are typically small in volume being incapable to 
host rare high-sigma peak objects. 

The requirement for large volumes and significant reso- 
lution to follow galaxy formation (gas inflows into relatively 
small scales) has made numerical studies of the growth of 
the first quasars extremely challenging. A number of ap- 
proaches have resorted to 'constrained' simulations. In a pi- 
oneering study, (Li et al. 2007) high-resolution merger sim- 
ulations with subgrid models for star formation and growth 
of black holes have been used to study the formation of 
z ~ 6 quasars. This work extracted merger trees from large, 
coarse dark matter only simulations and identified the most 
massive halo candidate at z = 6. To simulate at high resolu- 
tion the merger trees were populated with isolated galaxies 
which undergo the corresponding mergers and are then en- 
dowed with gas and models for star formation and black hole 
growth. 

This approach qualitatively reproduces the properties 
of the quasar SDSS J1148+5251 and its host at z = 6.42. 
Fully cosmological re-simulations of selected halos from the 
Millenium run were also carried out in a more recent study 
by Sijacki, Springel, & Haehnelt (2009). This approach also 
finds that the most massive halos are consistent with being 
the first quasar hosts. 

In order to study directly how and where the 
first quasars assemble, without a pre-imposed choice of 
halo, large cosmological volumes (~ Gpc 3 ) for capturing 
rare high-sigma peaks are required and sufficiently high- 
resolution in order to resolve kpc scales. High resolution 
is also necessary when including subgrid models for star- 
formation, black hole accretion and related feedback pro- 
cesses. Keeping these constraints in mind we have run the 
largest (currently feasible) hydrodynamic simulation of its 
kind, MassiveBlack , which includes gravity, hydrodynam- 
ics and subgrid models for star-formation, black hole growth 
and associated feedback processes. It was run in a cosmo- 
logical volume of L box = 533/i _1 Mpc with 2 x 3200 3 « 65 
billion particles (dark matter + gas) and a uniform gravita- 
tional softening of e = 5/i _1 kpc, with the code P-GADGET 



^box Afpart m DM m gas 6 

(fc^Mpc) (h^Mo) (H^Mq) (/i^kpc) 

533.333 2 x 3200 3 2.78xl0 8 5.65xl0 7 5 



Table 1. Basic simulation parameters for the MB simulation. 
The columns list the size of the simulation box, Lb ox > the number 
of dark matter particles used in the simulation, A^ par t, the mass 
of a single dark matter particle, m DM , the initial mass of a gas 
particle, m gas , and the gravitational softening length, e. All length 
scales are in comoving units. 

which has been extensively modified from the public code 
GADGET2 (Springel 2005) to run optimally on a large num- 
ber of multicore processors. The simulations were carried 
out on 10 5 cores of Kraken at NICS 1 . MassiveBlack has the 
same mass and force resolution (as well as similar subgrid 
models for star formation, black hole growth and associated 
feedback processes) as that of the resimulated halo of Li et 
al. (2007), however it does not rely on an imposed merger 
scenario to produce luminous quasars at z ~ 6 and tracks 
the assembly of galaxies in a more self-consistent manner. 

The high redshift observations of galaxies, quasars and 
their hosts can be used to test standard models of galaxy 
formation in previously unexplored regimes. Indeed Mas- 
siveBlack has been instrumental in reproducing a number 
of observational properties of high-z quasars, e.g. the forma- 
tion and abundances of Sloan Digital Sky Survey (SDSS) 
type black holes of mass M BH ~ 1O 9 M0 at z — 6 (Di Mat- 
teo et al. 2011) and also statistical properties, such as their 
luminosity function and high-redshift clustering (DeGraf et 
al. 2011). In this work we study the formation and proper- 
ties of galaxies hosting z ~ 6 quasars and as a sanity check 
we see how well observed global properties of galaxies com- 
pare with those in MassiveBlack . Our paper is organised as 
follows. We start with by describing the MassiveBlack sim- 
ulation in section 2. In section 2.1 we compare the GSMF 
in MassiveBlack with observations and earlier simulations. 
We next identify a sample of potential quasar host galax- 
ies at z r»j 6 and look at their formation and growth in 
sections 4.1 and 4.2. In section 4.3 we compare properties 
of these galaxies in MassiveBlack with recent observations. 
We present our conclusions in section 5. 



2 METHODS: THE MASSIVE BLACK 
SIMULATION 

In this section we describe a large hydrodynamic simula- 
tion, MassiveBlack , which we have run to study the high 
redshift Universe. We have used P-GADGET, a significantly 
upgraded version of GADGET3 (see (Springel 2005) for an 
earlier version) which we are developing for use at upcoming 
Petascale supercomputer facilities. MassiveBlack is a cosmo- 
logical simulation of a ACDM cosmology. 

It is worth pointing out that the number density of 
z rsj 6 quasars is extremely small, n ~ a few (Gpc) -3 , 
and that they are hosted by rare massive halos. There- 
fore in order to simulate and resolve these objects one 

1 http://www.nics.tennessee.edu 



Quasar Hosts at z ~ 6 3 



needs a large cosmological volume as well as high mass 
and force resolution. MassiveBlack was run with N pavt = 
2 x 3200 3 = 65.5 billion particles in a comoving volume of 
side Lbox = 533/i -1 Mpc and a comoving gravitational soft- 
ening length of e = 5/i _1 kpc, (see table 1 for more details). 

The initial conditions were generated with the Eisen- 
stein and Hu power ( spectrum at z — 159 and the simula- 
tion was evolved to z — 4.75. The cosmological parameters 
used were: the amplitude of mass fluctuations, as = 0.8, 
spectral index, n s = 0.96, cosmological constant parameter 
Ha = 0.74, mass density parameter Q m = 0.26 and baryon 
density parameter = 0.044. 

Along with gravity and smoothed particle hydrodynam- 
ics (SPH) P-GADGET incorporates a multi-phase ISM model 
with star formation (Springel & Hernquist 2003) and black 
hole accretion and feedback (Di Matteo, Springel, & Hern- 
quist 2005; Springel, Di Matteo, & Hernquist 2005). 

2.1 Subgrid model for black hole accretion and 
feedback 

Black holes are modelled as collisionless sink particles within 
newly collapsing halos in our simulation, which are identified 
by a friends-of-friends (FOF) (Davis et al. 1985) halofmder 
called on the fly at regular time intervals. A seed black hole 
of mass M S eed = 5 x 10 5 Mq is inserted into a halo with 
mass Mhaio > 5 x lO lo /i _1 M if it does not already con- 
tain a black hole. The seeding recipe is chosen to match 
the expected formation of SMBHs by gas directly collaps- 
ing to black holes with M BH ~ M see d (Bromm & Loeb 2003; 
Begelman, Volonteri, & Rees 2006) or by massive primordial 
stars collapsing into ~ 1O 2 M mass black holes (Bromm & 
Larson 2004; Yoshida et al. 2006) at z ~ 30. Once seeded, 
black holes grow by accreting surrounding gas or by merg- 
ing with other black holes. Gas is accreted with an accretion 

rate M BH = T° -%75 (Hoyle & Lyttleton 1939; Bondi & 

Hoyle 1944; Bondi 1952), where p is the local gas density, c s 
is the local sound speed and v is the relative velocity of the 
black hole and the surrounding gas. We limit the accretion 
rate to mildly super-Eddington consistent with Begelman, 
Volonteri, & Rees (2006); Volonteri & Rees (2006) to pre- 
vent artificially high values. The black hole radiates with a 
bolometric luminosity which is proportional to the accretion 
rate, Lboi = tiMbhc 2 (Shakura & Sunyaev 1973), where 77 is 
the radiative efficiency and its standard value of 0.1 is kept 
throughout, and c is the speed of light. Some of the liberated 
energy is expected to couple thermally to the surrounding 
gas. In the simulation 5% of the radiated energy does this. 
This energy is deposited isotropically on gas particles that 
are within the black hole kernel (32 nearest neighbours) and 
acts as a form of feedback energy (Di Matteo, Springel, & 
Hernquist 2005). The value of 5% is the only parameter in 
the model and was set using galaxy merger simulations (Di 
Matteo, Springel, & Hernquist 2005) to match the normali- 
sation in the observed M BH — a relation. Black holes also 
grow by merging with other black holes once they come 
within the spatial resolution with a relative velocity below 
the local gas sound speed. 

This model for the growth of black holes has been devel- 
oped by Di Matteo, Springel, & Hernquist (2005); Springel, 
Di Matteo, & Hernquist (2005). It has been implemented 



and studied extensively in cosmological simulations (Sijacki 
et al. 2007; Li et al. 2007; Colberg & Di Matteo 2008; Di 
Matteo et al. 2008; Croft et al. 2009; Booth & Schaye 2009; 
Sijacki, Springel, & Haehnelt 2009; Degraf, Di Matteo, & 
Springel 2010, 2011; Degraf et al. 2011; Chatterjee et al. 
2011) successfully reproducing basic properties of black hole 
growth, the observed M BH — a relation and the black hole 
mass function (Di Matteo et al. 2008), the quasar luminosity 
function (Degraf, Di Matteo, & Springel 2010) as well as the 
clustering of quasars (Degraf, Di Matteo, & Springel 2011). 

We identify galaxies with the group-finder code SUB- 
FIND (Springel et al. 2001). SUBFIND identifies bound 
clumps within a FOF halo and properties of galaxies such 
as position, mass (dark matter, gas, stars and black holes), 
star formation rate (SFR) among others are stored. 

We use a relational database management system devel- 
oped by Lopez et al. (2011) specifically for this simulation to 
track the history of black hole properties (e.g. mass, accre- 
tion rate, position, local gas density, sound speed, velocity, 
and black hole velocity relative to local gas) which are saved 
for each black hole at every timestep. For a complete sum- 
mary of the database format and its efficiency, the reader is 
referred to Lopez et al. (2011). 



3 RESULTS 

3.1 The Galaxy Stellar Mass Function 

We start by looking at the global properties of galaxies at 
high redshift in the MassiveBlack simulation. We use the 
halo catalogues generated by SUBFIND and pick all halos 
with more than 40 star particles to construct a galaxy stellar 
mass function (GSMF) from z = 7 to z = 5. This translates 
to a minimum stellar mass of galaxies of M* = 1O 9-2 M . 
We compare our GSMF with recent observed mass func- 
tions (Gonzalez et al. 2011) which were derived from the lu- 
minosity function at these redshifts (Bouwens et al. 2011a). 
We also compare our GSMF with those from recent hydro- 
dynamic simulations at these redshifts by Jaacks, Choi, & 
Nagamine (2011). We note that Jaacks, Choi, & Nagamine 
(2011) used multiple runs with varying resolution to achieve 
a larger dynamic range in resolved galaxies, and were run 
with the same SPH code and cosmological parameters as the 
MassiveBlack run in this study. Feedback due to AGNs was 
however not included by Jaacks, Choi, & Nagamine (2011). 
Two examples of the Jaack et al. runs are the N600L100 
(Lbox = 100/i _1 Mpc, iVpart = 2 x 600 3 ) of Jaacks, Choi, 
& Nagamine (2011) which has exactly the same mass res- 
olution as our MassiveBlack simulation and their N400L10 
run (Lbox = 10/i _1 Mpc, N part = 2 x 400 3 ) which can probe 
galaxies masses down to M* = 10 6 ' 8 M©. 

Our results (solid line) are shown in figure 1 for z = 5, 6 
and 7 (left to right). The data points are from observations 
at these redshifts (Gonzalez et al. 2011). The dashed lines 
are from Jaacks, Choi, & Nagamine (2011) and were only 
computed down to z = 6. At z = 5 our results are in good 
agreement with observations in the mass range that we are 
able to resolve. 

At z = 6 we slightly underpredict the GSMF for 
M* > 10 9 ' 5 M©, but a main difference with Jaacks, Choi, 
& Nagamine (2011) is that their GSMF curve drops off 
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Figure 1. The evolution of the GSMF from z = 5 to z = 7 (left to right) in the MB simulation (solid line), data points are from 
observations reported by Gonzalez et al. (2011). The dashed curve is from hydrodynamic simulations of high redshift galaxies by Jaacks, 
Choi, & Nagamine (2011). 



completely for M* > 1O 1U5 M . This may be due to the 
smaller volume (~ 150 times smaller than MB) that they 
have used. However they are consistent with our estimates 
for M* < 1O 1O M . 

At z — 7 both the theoretical GSMFs underpredict that 
of the observations in the entire mass range. The source of 
this discrepancy is at present unclear. This may be due to 
the relation between mass and light used by (Gonzalez et 
al. 2011) which was calibrated at z ~ 4 and is consistent 
with the sample at z ~ 5. At z ~ 6 and z ~ 7 this cor- 
relation is less obvious even though it is consistent, in zero 
point, with the even smaller sample at these higher redshifts. 
Our GSMF are evolving but the observations are consistent 
(within error-bars) with no evolution in the GSMF from 
z = 5 to z = 7. We will explore these issues as well as the 
luminosity function of high redshift galaxies, which is the 
direct observable in a future work. 



4 PROPERTIES OF HOST GALAXIES OF 
Z ~ 6 QUASARS 

We now focus on individual properties of host galaxies of z ~ 
6 quasars in this section. We look at how these galaxies and 
the black holes at their center were assembled and compare 
our results with their observed properties in the latter part 
of this section. 



4.1 Star formation and black hole growth 

Our simulation allows us to follow the growth of supermas- 
sive black holes and their host galaxies up to z ~ 5. In 
Figure 2, we plot the accretion rate of the black hole (black) 
and the SFR (red) of its host with redshift for a sample of 
eight luminous quasars selected such that their luminosity 
is consistent with the magnitude limit for quasars in the 
Sloan survey. This is the blue line, which is bounded by an 



z-band magnitude limit of m» < 20.2 for z > 3 (Shen et al. 
2009) and converted to a bolometric luminosity (and then 
an accretion rate) using the SED of Hopkins, Richards, & 
Hernquist (2007). This SDSS flux cut roughly corresponds 
to BH accretion rates of about 10 solar masses per year). 

In figure 2 we find that the star formation rate in the 
galaxy is strongly correlated with the growth of the central 
black hole. The central black holes grow rapidly through a 
period of sustained Eddington accretion ( typically between 
8 < z < 6) and are continuously fed by streams of high 
density gas (Di Matteo et al. 2011). The mass of the host 
halos in this sample grow from Mhaio ~ 10 11 " 6 at z = 8 to 
Mhaio ~ 10 12-4 at z = 6 . Prior to its peak accretion phase 
the black hole grows more rapidly that the stellar mass in its 
host. Star formation is regulated by feedback from the black 
hole, and typically star formation is suppressed just prior to 
the peak accretion phase of the black hole. This is a typical 
feature of this model and has been seen in many previous 
works (Di Matteo, Springel, & Hernquist 2005; Li et al. 2007; 
Sijacki, Springel, & Haehnelt 2009). Once the black hole ac- 
cretion becomes feedback dominated it deposits enough en- 
ergy in its vicinity and shuts off further growth by expelling 
gas in its surrounding star-forming region. At its peak, the 
SFR of the host galaxy is extremely high, C(10 3 )M o /yr , 
and is consistent with observations (Wang et al. 2010, 2011) 
for the SFR of quasar hosts at these redshifts. We make a 
direct comparison with observations in a later section. Most 
of the SMBHs fall within the Sloan detection limit when 
they grow through Eddington limited accretion and attain 
M BH = 1O 9 M0 during this phase. The one exception is the 
first object which peaks early on at z = 7.5 then undergoes 
a merger (as will become evident in figure 3) at z = 5.5 to 
become a M BH = 1O 9 M0 black hole. 

We look at the evolution of the environment around 
three typical host galaxies (the first three objects in the top 
row of figure 2) from z = 7.5 to z = 5.0 in figures 3-5. In the 
top row of these figures we plot the gas distribution, color 
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Figure 2. SFR(red) and black hole accretion rate(black) in M /yr. The SDSS flux limit of m, < 20.2 for z > 3 (Shen et al. 2009) for 
the quasar sample is shown by a blue line. 



coded by temperature. The middle row is the gas distribu- 
tion color coded by the SFR and the bottom panel is the 
distribution of stars. The large blue circle denotes the virial 
radius of the halo and the smaller circles are black holes with 
the radius proportional to the mass of the black hole. 

The SMBH in these examples have different growth his- 
tories, e.g. the first SMBH becomes feedback dominated at 
z = 7.5 then grows again due to a merger at z = 5.5 seen 
in figure 3. Apart from this SMBH the others grow to a 
mass of 1O 9 M by Eddington limited accretion by z ~ 6. 
As seen in figures 4-5 these halos are continuously fed by 
cold streams down to z ~ 6, prior to the feedback domi- 
nated phase. At this redshift feedback from the black hole 
starts to inject energy but is not able to disrupt the stream. 
Star formation is still sustained at SFR ~ lO 3 M /yr further 
down to z rsj 5.5. By z = 5.5 feedback from the black hole 
starts to destroy the inner structure of the cold streams, in- 
hibiting further growth; the SFR drops down by an order 
of magnitude of its peak value. There are still pockets of 
cold star- forming gas around the black hole, however the 
depleted gas due to star formation is no longer replenished 
by the cold streams. By z — 5 the inner part of the smooth 
cold stream is mostly destroyed, the gas in the central region 
is heated to T ~ 10 7 K and pushed out and dense subhalos 
which have been self-shielded from feedback contribute to 
the growth of the host galaxy and the central black hole, 
though at a much reduced rate. 

In the middle row of figures 3-5 we again look at the 
distribution of gas around the same object but now color 
coded by the SFR. The redshift of each panel is the same as 
in the top row. We find that the central region of the galaxy 
is forming most of the stars. Star formation also occurs in 
dense clumps of gas located on cold filamentary streams 



though at a much reduced rate. The central region has a 
sustained period of star formation down to z = 5.5 for most 
objects. During this period the black hole has grown at the 
Eddington rate through smooth accretion of cold gas. From 
this inspection it appears that most star formation also oc- 
curs in the most centrally-located gas that also feeds the 
black hole which, as demonstrated previously (Di Matteo et 
al. 2011) is cold-flow fed. 

At z = 5 the star formation has dropped drastically in 
the central galaxy due to feedback from the central black 
hole. It appears that feedback eventually destroys the cold 
streams inside the halos; star formation still persists though 
at a reduced level in dense clumps surrounding the black 
hole. In the bottom row of figures 3-5 we look at the distri- 
bution of stars. As expected the locations of stars are prefer- 
entially found near star- forming gas particles (middle row). 
We are unable to accurately predict the morphologies of the 
galaxies due to the limited resolution of the simulation. 



4.2 Growth of Host Galaxies Through Cold 
Streams 

Here we investigate the origin of the star forming gas in the 
quasar host galaxies. We will in particular test whether most 
of this star forming gas is indeed entering the halo via cold 
streams rather than cooling from a shock heated phase. We 
wish to point out that we do not characterise the redshift 
of accretion of these star- forming gas particles, but rather 
look at the origin of all the star-forming gas during the host 
galaxy's peak star- forming activity. 

As an example, we look at the object in figure 5. We 
examine how the star-forming gas in the halo at its peak 
star-forming activity, i.e. at z = 5.75, was accreted. We trace 
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Figure 3. An example of the growth of a typical z ~ 6 quasar host galaxy. This host corresponds to the first object in the top panel of 
figure 2. The top row visualizes the gas distribution color coded by temperature across six redshifts. The projected density ranges from 
~ 10 5 to ~ 1O 8,5 M0 (kpc/h) -2 . In the middle row we visualize the gas distribution but now color coded by the SFR. The bottom row 
shows the distribution of stars. The large blue circle is the virial radius of the halo. The smaller circles are black holes with the radius 
proportional to their mass. 
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Figure 4. Same as in figure 3 but for the second object, top row, in figure 2. 



the temperature histories of these star-forming gas particles 
back to z = 7.5 and plot them as a function of (physical) 
separation from the SMBH, in figure 6. To make compari- 
son with observations easier, in the discussion that follows 
all length scales are quoted in physical units. The temper- 
atures represent the effective temperature for star-forming 
gas particles which are in the two-phase medium (Springel & 
Hernquist 2003) and the temperature for non star-forming 
gas which are in the single phase medium. 



In figure 6 the filled black circles denote star-forming 
gas particles and the open blue squares denotes the gas par- 
ticles which have zero SFR. The horizontal dot-dashed line 
is the virial temperature, the vertical solid line is the virial 
radius and the vertical dashed line is the gravitational soft- 
ening length. Prior to z = 5.75, one can identify four distinct 
regimes in the T — r plots of the gas particles which end up 
forming stars in the halo: (1) a cold stream of non- star form- 
ing gas with temperatures of ~ 10 4 K beyond r ~ lOkpc, (2) 
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Figure 5. Same as in figure 3 but for the third object, top row, in figure 2. 



clumps of star- forming gas outside r ~ lOkpc with tempera- 
tures 10 4 < T < 10 6 K, (3) hot non-star forming gas outside 
r ~ lOkpc with temperatures 10 5 5 < T < 10 7 K and finally 
(4) dense star-forming gas within r ~ lOkpc with tempera- 
tures 10 4 < T < 10 6 K. Star formation only occurs in dense 
environments, i.e. at the centre of the halo and clumps lo- 
cated in filaments, this is also seen in figures 3-5. A fur- 
ther investigation reveals that the cold non star-forming gas 
(with T ~ 10 4 K) is also located in filaments whereas the 
hot non star-forming gas is diffuse and spread out across 
the halo. 

As can be seen in figure 6 most of the gas that is 
star- forming at z = 5.75 does not come from the diffuse 
hot medium. The major mode of gas accretion for the host 
galaxy is from gas in filamentary streams. The streams have 
two components, dense star forming clumps and cold non 
star-forming gas that penetrates deep into the halo well 
within the 10 kpc region of the central black hole. At this 
point gas is dense enough to form stars in the central galaxy 
and also fuel the growth of the black hole. The mass of 
gas in star forming clumps in the stream is small compared 
to the cold non- star forming gas in the stream. Feedback 
from supernovae and the black hole heats up the gas at the 
centre, so that the temperature systematically increases to 
T ~ 10 6 K with decreasing separation. This mode of accre- 
tion continues down to z ~ 5.75. 

Eventually by z = 5 the black hole has injected enough 
energy into the surrounding medium to destroy the cold 
stream supplying gas to the central galaxy. The black hole's 
growth has become self-regulated. The black hole has also 
regulated the growth of its host. We find little star- formation 
within 3 kpc of the black hole, but a residual amount of dense 
gas is still clumped in the region with 3 < r < 10 kpc, so 
that star-formation still persists in this part of the central 
galaxy. 

A sustained period of cold accretion is largely respon- 
sible for the high SFR for most objects. The exception is 



the first object in figure 2 where a merger at z = 5.5 is re- 
sponsible for increasing the SFR. As seen in figures 3-5 star 
formation occurs only in dense regions, i.e. mostly in the 
central halo and at a reduced level in dense clumps within 
filaments. At its peak, 95% of star formation occurs within 
the 3 kpc region of the central galaxy. The numbers quoted 
for this example are representative of those for the full sam- 
ple, fluctuating by only a few percent for the other galaxies. 

We track the entire temperature history of the star- 
forming gas particles of the host galaxy at its peak star- 
forming phase, i.e. at z = 5.75. We define Tmax as the max- 
imum temperature that a gas particles attains prior to its 
first star-forming phase, i.e. prior to it first being in a two- 
phase medium (Keres et al. 2005, 2009). A distribution of 
Tmax will then indicate whether these gas particles were 
accreted through the hot or cold mode (Keres et al. 2005, 
2009). This is shown in the bottom right panel of figure 6. As 
can be seen the majority (> 85%) of the star- forming gas is 
accreted through the cold mode, the remaining gas (< 15%) 
comes from gas which has been heated while accreting onto 
the halo and then cools to form stars. These hot gas par- 
ticles are the same as the diffuse hot non star-forming gas 
seen in the first four panels of figure 6. The numbers pre- 
sented here represent a lower bound for cold mode accretion 
for star-forming gas, since a larger fraction of star-forming 
gas particles prior to z = 5.75 should have accreted onto 
the halo through the cold mode. We find that if we consid- 
ered all the gas, not distinguishing between star- forming and 
non star-forming gas, then ~ 70% of the gas at z = 5.75 is 
accreted cold and ~ 30% is accreted through the hot mode. 

The analysis presented here indicates that the quasar 
host galaxies are forming stars very efficiently. The massive 
star formation is sustained through a supply of gas mostly 
from cold streams, that are able to penetrate the halo and 
reach the center without being heated to the virial tem- 
perature of the halo consistent with earlier work at lower 
redshift (Dekel et al. 2009). A smaller fraction of gas < 15% 
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Figure 6. Temperature histories of the star-forming gas of the central halo at z = 5.75 as a function of separation from the central 
SMBH. All length scales in this figure are in physical units. This particular halo is the same as in figure 3 and has a peak SFR at 
z ~ 5.75. Star forming gas particles in the central halo are tracked backwards with respect to the reference redshift of z = 5.75 when this 
halo was at its peak star- forming activity. The horizontal dot-dashed line denotes the virial temperature (Kelvin) of the halo, the solid 
horizontal line is the virial radius of the halo and the dashed horizontal line is the gravitational softening length. The filled black circles 
denote gas particles with non-zero SFR and the open blue squares are gas particles with zero SFR. The bottom right panel indicates the 
distribution of T macc for star- forming gas particles at z = 5.75 



does get heated to ~ T v - ir to finally cool and form stars. At 
its peak ~ 95% of the star- format ion is occurring within a 3 
kpc region of the central halo, consistent with observations 
(Walter et al. 2004, 2009; Wang et al. 2010). 

4.3 Comparison with Observations 

In this section, we compare properties of quasar hosts in the 
MassiveBlack simulation with recent observations (Wang 
et al. 2010). These observations specifically constrain ob- 
servables for the host galaxies such as SFR, molecular gas 
M mol = M [H 2 + He], and the M BH - a relation. Addition- 
ally we look at the M BH — M* relation for the host galaxies in 
the MassiveBlack simulation and compare it with the local 
relation in observations. 

4-3.1 Star Formation Rates and Cold Gas 

The reprocessed emission in the far-infrared (FIR) from star 
formation-heated dust is used to provide an estimate of the 



star formation from observations of z ~ 6 quasar host galax- 
ies (Wang et al. 2008). The AGN contribution is removed 
and the remaining FIR luminosities, Lfir, are then con- 
verted to a SFR (Wang et al. 2010, 2011), using 




The sample of quasar host galaxies in the observa- 
tions of Wang et al. (2010, 2011) fall in the redshift range 
5.78 < z < 6.43. We compile the SFR of the quasar hosts in 
our simulation (the ones in Figure 2) and compare it with the 
observations in the left panel of Figure 7. We find that most 
of these host galaxies have a peak SFR ~ lO 3 M0.yr _1 in 
the redshift range 5.5 < z < 6.5, comparable to the observed 
SFR of quasar host galaxies. However our model seems to 
lack objects close to the observed SFR of 15OOM0yr -1 for 
some host galaxies. It is of course possible that our star for- 
mation model is somewhat simplistic, particularly at these 
redshifts, since we do not model star formation from molec- 
ular gas (e.g., Krumholz & Gnedin (2011)). Further investi- 




Figure 7. Le/£: Evolution of SFR for the hosts of the luminous quasars. Right: Evolution of M cold (solid lines) and M mol (dot-dashed) 
of the hosts of the most massive black holes. Data points in both panels are from observations (Wang et al. 2010, 2011). Data points on 
the right panel are estimates for M mol . 



gation of the star formation modelling is beyond the scope 
of this paper so instead we prefer to compare directly to as- 
sociated measurements of the cold and molecular gas, which 
is what the observations can constrain. 

In the right panel of figure 7 we plot the cold gas mass, 
M cold , of the host. We find that the the evolution of M cold 
broadly follows the trend seen in the SFR. Observations on 
the other hand probe the molecular gas mass, M mol of these 
hosts which are not directly modeled in our simulations. 

The size of the cold molecular gas reservoir which fuels 
star formation in these observed quasar host galaxies has 
been estimated through redshifted CO emission, (Bertoldi 
et al. 2003; Walter et al. 2003; Carilli et al. 2007; Wang 
et al. 2010). These studies indicate molecular gas masses 
of > 1O 1O M in these objects. Since our simulations do not 
model molecular gas instead estimate the amount of molecu- 
lar gas by using the SFR as a proxy for M mol , i.e. by convert- 
ing the SFR to Lfir, eq. 1 and using the relation between 
the CO (1-0) line luminosity, ^co(i-o)5 an< ^ ^fir for local 
star- forming systems such as local starburst spiral galaxies, 
Ultra Luminous Infra Red Galaxies (ULIRGs) and high-z 
submillimeter galaxies (SMGs) (Solomon & Vanden Bout 
2005). 

log (Lfir) = 1-7 x log (Lco(i-o)) - 5 (2) 

Lqo(i-o) can then be converted to a cold molec- 
ular gas mass M co id = aL / CO ( 1 _ ^ with a = 
O.8M (K km s" 1 pc 2 ) -1 (Wang et al. 2010). Using these 
relations we look at the evolution of molecular gas for quasar 
host galaxies in the MassiveBlack simulation and compare 
them with observations (Wang et al. 2010, 2011) in the right 
panel of figure 7 shown as dot-dashed lines. Here again we 
find that we are able to reproduce the amount of molecular 
gas in these galaxies at redshifts 5.5 < z < 6.5. The esti- 



mates are in better agreement than the observed SFR. Wang 
et al. (2010) also find that z ~ 6 quasar hosts lie above the 
local Lfir-Lco(i-o) relation and attribute it to unknown 
AGN contributions to Lfir and different dust temperatures 
and CO line ratios, even though they attempt to correct for 
AGNs. 

The ratio between M mol and M cold is typically 10%. 
This is much higher than 1% seen in observations of local 
galaxies (Obreschkow & Rawlings 2009). This indicates that 
z ~ 6 quasar host galaxies are forming stars mores efficiently 
than local galaxies. 

4.3.2 The M BH - a and the M BH - M* relation 

We now look at the evolution of M BH — a and the M BH — M* 
relation between the central black hole and its host galaxy, 
where a is the velocity dispersion of the bulge, for which 
we use the velocity dispersion of stars within the half-mass 
radius as a proxy. Given that the black hole grows faster than 
its host, as seen in figure 2 we expect to see these relations 
to be steeper than the local relations, e.g. Tremaine et al. 
(2002); Haring & Rix (2004). 

The evolution of the M BH — a relation (solid line) is 
shown in the left panel of figure 8 for the MassiveBlack sam- 
ple. Comparison is made to the best fit local relation (dot- 
ted line) (Tremaine et al. 2002) and observations of z ~ 6 
quasar-host systems (triangles) (Wang et al. 2010). In this 
panel we show the evolution up to z = 5.5 (solid lines) in 
order to better compare with observations at these redshifts 
and further down to z = 4.75 (dot-dashed lines). Since the 
inclination angle is degenerate with the CO linewidths, an 
average inclination angle of 0i nc = 40° was assumed (Wang 
et al. (2010), filled triangles); the open triangles relax this 
constraint. As expected we see that the black holes grow 
more rapidly than their host galaxies and eventually end up 
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Figure 8. Evolution of Mbh — & (left) for quasar-host systems in the MassiveBlack simulation up to z = 5.5 (solid line) and further 
down to z = 4.75 (dot-dashed line). The dotted line is the best fit local relation of Tremaine et al. (2002). The filled triangles are data 
from observations z ~ 6 quasars and host galaxies (Wang et al. 2010) assuming an average inclination angle 9- luc = 40° and open triangles 
are the same observations without any assumption for the inclination angles. Evolution of the Mbh — M* relation (right) for the same 
objects as in the left panel. Solid lines follow the evolution to z = 5.5 and dot-dashed lines continue it down to z = 4.75. Comparison is 
made with the the best fit local relation of Haring & Rix (2004). 



above the local relation of (Tremaine et al. 2002). By z = 5.5 
the M BH — a relation compares reasonably well with observa- 
tions which assumes an inclination angle. However observa- 
tions indicate that there is an object with a ~ 120km/s and 
M BH ~ 1O 9-7 M0 even when the assumption of inclination 
angle is considered. We do not find such an object within 
our sample and this discrepancy may have to do with the 
uncertainties in the inclination angle. 

Next we look at the evolution of the M BH — M* rela- 
tion in the MassiveBlack simulation (solid lines) and com- 
pare with the best fit local relation of Haring & Rix (2004) 
(dotted line) in the right panel of figure 8. Again consistent 
with previous properties such as the SFR-Mbh relation and 
the M BH — a relation we find that the black hole is assem- 
bled more rapidly than the central host galaxy so that the 
M BH — M* relation is well above the local relation of Haring 
& Rix (2004). It would be interesting with future simula- 
tions to see what mechanisms eventually drive them to the 
local relation. 

The steep growth of M BH which outpaces the stellar 
component of the host galaxy is common for all these sys- 
tems and is again consistent with the fact that the accretion 
rate of the black hole as a function of redshift has a steeper 
slope than the SFR. The cyan and red lines in figure 8 corre- 
spond to the first two objects of in the top row of figure 2. A 
closer inspection shows that these objects undergo relatively 
dry mergers with subhalos in the feedback dominated phase. 
The subhalo has a larger stellar component compared to gas 
and therefore during this merger both a and M* grow more 
rapidly than M BH . This can be seen as a relative flatten- 
ing of the M BH — a and M BH — M* relation in figure 8. It 
is plausible that dry mergers are responsible for eventually 



bringing the high redshift M BH — a and M BH — M* rela- 
tions onto the well constrained local relation (Tremaine et 
al. 2002; Haring & Rix 2004) 



5 CONCLUSIONS 

MassiveBlack has been successful in reproducing a number 
of important properties of quasars at high redshift within 
observational constraints: (i) the formation and abundances 
of Sloan-type black holes of mass M BH - 10 9 M o at z = 

6 (Di Matteo et al. 2011), (ii) statistical properties such 
as the luminosity function of quasars and its evolution as 
well as their high-redshift clustering (DeGraf et al. 2011). 
Importantly MassiveBlack has indicated that cold flows can 
sustain the growth of black holes at Eddington rates so that 
they attain masses of ~ 1O 9 M0 in a short time span (Di 
Matteo et al. 2011). 

In this paper we have looked at the formation of galaxies 
hosting these z ~ 6 quasars with the MassiveBlack simula- 
tion. We are able to reproduce a number of observational 
properties of these galaxies. We summarize our findings be- 
low: 

• We find good agreement between the theoretical GSMF 
produced by MassiveBlack and observations at z — 5 and 
z = 6. 

• At z ~ 6, M BH rsj 10 9 M© black holes are already in 
place. The growth of the black hole and its host galaxy are 
strongly correlated. The black hole regulates the growth of 
its host galaxy. 

• Black holes grow faster than the host galaxy and this 
is reflected in the deviations seen in our study from the 
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local M BH — a relation and M BH — M* relation; our results 
are however consistent with observational findings for these 
deviations at z ~ 6. 

• The cold streams that feed the black hole are also re- 
sponsible for feeding the host galaxy to sustain the excessive 
starburst in them. The accretion of gas, which eventually 
form stars in the host, through the hot mode is subdomi- 
nant ( < 15%), the dominant fraction being accreted cold 
during their history. 

• Host galaxies grow in extreme environments. We are 
able to reproduce the high SFR of ~ lO 3 M yr _1 seen in 
observations in redshifts 5.5 < z < 6.5. Our simulation sug- 
gests that the observed galaxies are at the peak of their 
star- for mat ion activity at these epochs. 

• Our derived estimates of molecular gas in these galaxies 
are consistent with observations. Such large reservoirs of cold 
molecular gas, M cold > 1O 1O M are responsible for fuelling 
star-formation in these galaxies. 

• Using the SFR as a proxy for M mol we find that the 
ratio of molecular and cold gas is larger than seen for local 
galaxies. This again indicates the host galaxies are forming 
stars more efficiently than their local counterparts. 

• We find that most of the star formation occurs within a 
compact 3-6 kpc (physical) region around the black hole at 
the peak of its star-forming activity. These scales are again 
consistent with observations (Walter et al. 2004, 2009; Wang 
et al. 2010). 

• Given that we are able to reproduce observations be- 
tween 5.75 < z < 6.5 our simulations suggests that these 
quasars and their host galaxies are seen at the peak of their 
assembly. 

From the results presented here we expect that com- 
pared to the local Universe the M BH — a relation is very 
different in the high redshift Universe. This has been indi- 
cated in earlier work (Di Matteo et al. 2008), though the 
sample was too small to put strong constraints on the rela- 
tion at z — 5 and above. Given the large volume and high 
resolution (hence large sample of galaxies and black holes) 
MassiveBlack is well suited for predicting the evolution of 
the M BH — a and M BH — M* relations at z = 5 and above. 
It would also be interesting to look at the mechanisms re- 
sponsible for their evolution. We will address these issues in 
subsequent work. 
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